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Abstract 

Despite its fundamental and technological importance, a microscopic understand- 
ing of the crystallization process is still elusive. By computer simulations of the 
hard-sphere model we reveal the mechanism by which thermal fluctuations drive the 
transition from the supercooled liquid state to the crystal state. In particular we 
show that fluctuations in bond orientational order trigger the nucleation process, con- 
trary to the common belief that the transition is initiated by density fluctuations. 
Moreover, the analysis of bond orientational fluctuations shows that these not only 
act as seeds of the nucleation process, but also i) determine the particular polymorph 
which is to be nucleated from them and ii) at high density favour the formation of 
fivefold structures which can frustrate the formation of crystals. These results can 
shed new light on our understanding of the relationship between crystallization and 
vitrification. 
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The liquid-to-solid transition is characterized by the spontaneous breaking of both po- 
sitional and orientational symmetry, but how this happens microscopically is still a matter 
of debate jJ-lSl]. Most approaches, like classical nucleation theory (CNT) or density func- 
tional theories (DFT) j?, 8|, assume that the crystallization process is primarily controlled 
by positional ordering, with the liquid regarded as a spatially uniform background where 
nucleation can occur at any location with an equal probability. However experiments 



and simulations 



mm 

have recently started to point out deviations from the classical pic- 



ture of crystallization, suggesting that this process could be more complex than previously 
thought. 

We argue that for understanding the origin of such deviations it may be crucial to rec- 
ognize the role of thermally excited fluctuations in driving the transition from the liquid 
phase to the crystal phase. Fluctuation effects were first identified in globular proteins and 
colloidal systems close to a metastable critical point, where crystallization starts with the 
formation of amorphous high-density aggregates and is followed by the actual nucleation 



event occurring within these fluctuations 



16 



-|20|: the two-step nucleation scenario. These 



studies revealed that the coupling between critical concentration fluctuations and density 
ordering (crystallization) plays a key role in nucleation. Even for a single component liquid. 



experiments 



ll| and simulations [1^, [ij] have recently showed the importance of density 



fluctuations in the initial stage of crystallization, which leads to the formation of precursors. 
Since the two-step nucleation scenario looks valid far j2o| or even in absence jl^ of a critical 
point, it has been suggested that this scenario (in which density fluctuations foreshadow 
structural ordering) could indeed be a general nucleation mechanism. Independently from 



13 



2l| have pointed out the 



the aforementioned two-step scenario, recent simulation works 

importance of another type of fluctuations occurring in the supercooled liquid phase: spon- 
taneous critical-like fluctuations of bond orientational order 



22 



23| . While the density order 



parameter (and in general translational order) is a measure of the relative spacing between 
the neighbouring particles, bond orientational order expresses instead the relative orienta- 
tion of the (geometrical) bonds between a particle and its neighbouring particles. In both 
scenarios, thermal fluctuations promote the formation of crystal precursors, i.e. preordered 
regions which trigger the nucleation process. However, since density and bond orientational 
ordering proceed simultaneously in the process of crystal nucleation, it has remained elusive 
how these order parameters are coupled, and whether any of the two plays a primary role. 
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In the present work we will investigate precursors in models of colloidal systems in or- 
der to elucidate the microscopic mechanism of crystal formation. We use here the word 
precursor as a short term for denoting the region of the liquid's free energy basin where 
nucleation is more likely to occur. We will first rule out the possibility of a two-step process 
involving densification as the first step towards crystallization. We will show instead that 
the nucleation process proceeds with the crystalline structures emerging first at liquid-like 
densities, a process akin to what was reported by some studies of nucleation in molecular 
systems 



24 



25|. By examining the crystallization process in the two dimensional order- 



parameter space of density and orientational order, we will show that precursor regions are 
not characterized by locally denser regions, but by locally bond-oriented regions, and we 
will present a novel microscopic explanation of this mechanism. We will show that these 
precursor regions not only act as seeds of the nucleation process, but also determine the 
particular polymorph which is to be nucleated from them. This new concept implies that 
polymorphism is already a property of the metastable liquid state. 

It is interesting to note that regions of high bond orientational order have also been identi- 
fied as responsible for the highly heterogeneous dynamics in deeply supercooled liquids, and 



could be linked to a growing structural length at the origin of the glass transition [22|, |26 |. 
A study of the microscopic properties controlling the crystallization of the liquid is thus of 
utmost importance not only in elucidating the pathway to crystallization, but potentially 
also to explain how crystallization can be avoided. In this context, we will show that our 
two-order parameter description provides a thermodynamic justification of Frank's hypoth- 
esis [27] that icosahedral clusters of particles act as inhibitors of crystallization. 

In this Article we concentrate on the homogeneous nucleation process for the simplest 
nontrivial model of a liquid, hard spheres of diameter a, by means of computer simula- 
tions. This system is ideal for studying crystallization and has already provided tremendous 



contributions to our basic understanding of crystal nucleation 28 



30[. In Supplementary 



Information we extend the generality of our study by applying the same concepts to very 
different classes of materials, in particular systems governed by ultrasoft potentials (like 
polymeric materials) and tetrahedrally coordinated potentials (like water). 
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Results 



Let us start by introducing the order parameters used in this study. We will describe here 
their basic properties, while for the exact mathematical definition we refer to the Methods 
section. We will always adopt a microscopic approach, by studying local order parameters 
(defined at a particle level). Since the liquid-to-solid transition is characterized by both 
translational and orientational symmetry breaking, we wish to monitor both properties 
during the crystallization process. A good order parameter for translational order, which 
expresses the relative spacing between particles in the system, is of course the local density 
Pi. This is easily computed by means of Voronoi diagrams, which assign to each particle 
a local volume vi = 1/ pi. To describe orientational order, which expresses the relative 
orientation between the neighbours around each particle, we use the spherical harmonics 



analysis introduced by Steinhardt et al. [Sll- We thus define our bond orientational order 
parameter as qQ{i), which is a rotationally invariant scalar defined for each particle i. A 
closely related order parameter is Qq{i)., which is obtained by coarse-graining qQ{i) over its 
neighbours. The importance of Qq lies in the fact that it is a good order parameter to 
detect precursor regions, as we will show later. Finally, to address the question whether 
crystal nuclei emerge from dense precursors, we need an order parameter that distinguishes 
disordered configurations from crystal-like ones. We call this order parameter S (as for 
structure): for particle i it goes from a value close to in the liquid-state to a value close 
to 12 (as the number of neighbours in a close-packed structure) in the crystal state. 

Composition of crystals during nucleation and growth 

We begin by following 50 spontaneous crystallization events from the metastable state at 
reduced pressure [ipa^ = 17, where (3 = l/ksT and a is the hard-spheres diameter. Under 
these conditions nuclei form and dissolve repeatedly, until the appearance of a nucleus which 
grows over the critical size and eventually spans the whole system. For each configuration 



we identify crystal particles following the criteria pioneered by Frenkel and co-workers [28 1 
(see Methods), and identify individual clusters via a cluster algorithm. Figure [1^ shows the 
average number of particles with local bcc, hep or fee coordination within the crystal nuclei, 
as a function of their size. The vertical dashed line in Fig. [1^, which indicates the average 



4 



size of the critical nucleus (ric ^ 80) obtained from umbrella sampling simulations (see 
Supplementary Information), separates the nucleation and growth regime. Within clusters 
of size smaller than ric, (66 ± 1)% of the particles are in local fee coordination. This is 
markedly different from the ratio for random stacking of hep and fee hexagonal planes, 
rifcc/T^hcp ~ 1) which is predicted from the very small bulk free energy difference (around 0.1% 
of the thermal energy in favour of fee) between fee and hep phases 32|, |33| . This behaviour 



of hard spheres, also pointed out in earlier studies including both experiments [2|, |29| 



and simulations 



37 



34-: 



remains to our knowledge still unexplained and we will show in the 
following a mechanism which accounts for this unbalance. The inset of Fig. [T^ shows the 
average density of the crystalline particles as a function of the nucleus size. All crystalline 
phases form at an average number density of ~ 1.06cr~^, higher than the metastable liquid 
density of ~ 1.02o'~^. The presence of a jump is of course expected for the averaged order 
parameters (both p and q^) ai a. first-order phase transition. More surprisingly instead, the 
density at which the smallest crystals start forming is still very far from the bulk density of 
the stable crystal {ps — 1.136cr~^). Thus the nucleation of the solid phase happens under 
conditions very far from the bulk solid. As the crystal grows, both the densities of the fee 
and hep phases gradually increase, whereas bee particles are unable to pack efficiently, and 
hence do not contribute to the cluster growth. Here we note that a bulk bcc crystal is in fact 
mechanically unstable in hard spheres (meaning that a bulk bcc crystal will immediately 
transform into a mixture of fee and hep crystals). 

Now we turn to the order parameter profiles of crystal nuclei. Figure Hb shows the 
averaged radial profiles of p(r) for different sizes of the nucleus (indicated by the arrow). 
The density profiles gradually increase as the nucleus becomes bigger, but still do not reach 
the bulk values even for sizes much larger than the critical nucleus size. This is in stark 
contrast to the prediction of classical nucleation theory (CNT), according to which critical 
nuclei share the same thermodynamic properties of the bulk solid phase. Such deviations 



from CNT was predicted by non-classical approaches 



IJ4QI 
2q|, 



-|42|. Contrary to a two-step 

scenario, where densification foreshadows structuring [iJ, l20|, we find no such an indication, 
as shown in the inset of Fig. [TJ), where the density gap Ap between the nucleus and the 
liquid phase is displayed for different radii R/Rcviticax (normalized to the value of the critical 
radius). The density of the nucleus grows continuously from the liquid, with an almost linear 
relationship between Ap and the nucleus size R. 
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Figure [U: shows both the density radial profile p{r) and the profile of the structural order 
parameter S{r) for critical nuclei (n ~ 80). Both profiles are normalized as to be unity in 
the pure fee crystal, and zero in the bulk liquid phase. Going from the liquid phase (r = oo) 
to the centre of the nucleus (r = 0) we see that the nucleus first develops some structural 
order at liquid-like densities, and only later does the density increase as well. At the centre 
of the nucleus both the structural order parameter and the density are far from their bulk 
values, but density is lagging behind the development of structural order. The inset shows 
the (S',p)-map for nuclei of different sizes. The continuous line is the classical behaviour, 
while simulation points always fall in the region of structured precursors, and not locally 
denser precursors. We note that the gradual increase of structural order is rather similar 



to that reported in 



24l |. where the structural order profile grows both its height and range 



24| is derived from a one- 



simultaneously. It may be worth noting that the result in Ref. 
order-parameter DFT model, where a perfect decoupling of structural order from density is 
implicitly assumed. The introduction of a coupling between density and structural order in 
the same type of model leads instead to the saturation of both structural and density order 
at the first stages of nucleation 17|, |40| . This is an interesting point to be studied since, as 
described later, our results suggest indeed a weak coupling between the two types of order 
parameters. In relation to this, we also note that translational order in DFT is not the same 
as bond orientational order: the former is specific to solid-type fluctuations, but the latter 
can be linked to both liquid-type and solid-type fluctuations. 

In conclusion we have found no signs of the two-step process involving enrichment at 



constant size and then growth, contrary to some theoretical predictions IJ, l20|, |42|. We 



rather find that the density increase is foreshadowed by the prestucturing of the nucleus. 
This lagging of densification behind structuring is similar to the results of previous nucleation 
studies in Lennard- Jones systems js], l^,^, 40|, but with the difference that in these studies 
both density and structural order are already saturated to the equilibrium values when the 
nucleus size slightly exceeds the critical size, whereas not in our case (see Fig. [Tb and 
c). Moreover the prestructuring prior to densification has always been ascribed to the low 
compressibility of the liquid phase (see, e.g., jsj). In the next section we will show instead 
that density fluctuations in the liquid and crystal phase overlap to a large extent, and that 
the prestructuring of the nucleus is rather due to the development of orientational order, as 
the true first step towards crystallization. 



6 



Interplay between density and bond orientational order 



To explain the nucleation pattern unveiled in the previous section, we will address the 
question of how density (p) and orientational order (gg) are coupled. In Fig. we display 
{p,qe)-ma.ps for the metastable liquid (before the appearance of the critical nucleus) at differ- 
ent pressures. We average separately for particles identified as liquid (liquid branch, dashed 
line) and crystal (crystal branch, lines with symbols) (see also Fig. Sib in Supplementary 
Information). By comparing the relative position of the two branches in the {p,qQ)-map 
it is easy to spot the regions of stability of each phase: the stable branch lies below the 
metastable branch, having higher orientational order at fixed density (or conversely, the 
stable phase can reach the same degree of orientational order at lower packing). Let us 
start by examining the system at reduced pressure Ppa^ = 11. This pressure is just below 



the melting pressure, which is Ppa^ = 11.54 43|. As shown in Fig. the liquid branch is 



always located below the crystal branch, and it is thus the stable branch for all values of 
Qq and p. This result is of course the expected one, since we are before the melting line. 
What is surprising is that we are able to determine the relative position of the system with 
respect to the melting line by simply looking at its (p, q^) map, without resorting to free 
energy calculations. And again as expected, as we increase the pressure a crossover between 
the two branches appears, with the crystal branch gaining stability. For clarity we will focus 
on the curves at (3pa^ = 17, which is the same pressure at which we obtained our Fig. [H 
At low p and gg the liquid branch is the stable one. The crystal branch remains metastable 
until it reaches a plateau of constant p, where the crossover with the liquid branch occurs. 
The value of this plateau is p = 1.06a~^ which is exactly the average density of the onset of 
crystal formation which we determined in the inset of Fig. [1^. After this plateau the crystal 
phase thus becomes the stable phase. This means that in the metastable liquid, particles 
which reach (because of thermal fluctuations) values of Qq and p bigger than the crossover 
values are in local coordination shells that are transforming from liquid to crystal-like. The 
reason why this process occurs at constant density is clear if we consider the fact that these 
particles are already embedded in regions of high orientational order. This means that their 
neighbours are already highly ordered, and by means of small local rearrangements are able 
to attain the symmetry of the crystal (in practice crossing the threshold which we use to 
identify crystal particles). We show an example of such a microscopic process in the snap- 
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shots in Fig. |2^, where small local rearrangements (white arrows) cause a change of the 
coordination around the central particle from liquid-like (blue) to crystal- like (red), without 
changing the local density but by increasing significantly the orientational order. 

If we now follow the curve at higher gg and p a surprising result emerges: a second 
crossover between the crystal and liquid branches makes the liquid branch stable again. The 
density of this crossover is p = l.lOTcr"^, which corresponds to a volume packing of ~ 58% 
(which is also the conventional value which marks the beginning of the glassy state in hard 
spheres 30|). This second crossover tells us that at very high and p the crystal phase 
becomes unfavoured again. Note that these are purely static results, not affected by the 
underlying dynamics. By using bond orientational analysis (see Supplementary Informa- 
tion), the structures responsible for the stability of the liquid branch at high density are 
easily identified as particles embedded in icosahedral environments. Icosahedral particles 
belonging to the liquid branch can attain higher densities than the corresponding crystal 
structures, but due to their fivefold symmetry are not able to attain long range translational 
order. The second crossover in the (p, ge) map tells us thus that crystals have a stability 
window, which is limited at low densities by disordered configurations (larger configurational 
entropy of liquid particles), and at high density by clusters with icosahedral structure. We 
have thus shown that icosahedral particles act as inhibitor to crystallization, as was recently 



observed in both experiments 



44, 



45| and simulations j46|. This is consistent with a sce- 



nario that glass-forming ability is controlled by frustration against crystallization, or the 



presence of 



a liquid 



23 



ow 



26 



ree-energy local configurations incompatible with the crystal symmetry in 



47- 



We have seen that the crystallization process is driven by the development of orientational 
order, which explains the prestructuring of the nuclei at liquid-like densities. Precursor 
regions are thus easily identified by bond orientational order alone. The one disadvantage 
of ge is that it also reveals the signal from icosahedral environments of particles. To locate 



crystal precursors, an effective strategy is to spatially coarse-grain qq 13 



50| , thus enhancing 



the signal from coherent regions (crystal-like) and suppressing it in disordered or icosahedral- 
like regions. This is the order parameter called Qe? which grows continuously from the liquid 
branch to the crystal branch. In Fig. [2]d we plot, for the metastable liquid (prior to the 
appearance of the critical nuclei) at pressure Ppa^ = 17, a map in the (Qe, p) plane of the 
structural order parameter S. S{i) quantifies how many first-shell neighbours of particle i 
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have similar local environments: for a disordered liquid we expect 5* to be null, whereas for 
a bulk close-packed crystal to be 12, i.e. all neighbours share the same environment. As 
we can see from Fig. the structural order parameter grows 'continuously' from low Qq to 
high values. Contour lines are almost parallel to the p axis, meaning that density is only 
weakly coupled to the increase of crystalline structure. In other words, high density regions 
encompass all possible values of the while high Qq regions are always the most crystalline. 
So precursor regions are exclusively controlled by the coarse-grained orientational order 
parameter, and density fluctuations are not sufficient to promote crystallization. 

Polymorph selection 

Crystals repeatedly appear, grow and melt as represented by the fluctuations in the 
bond orientational order parameter Q^. Since crystal nuclei appear from regions of high 
bond orientational order, the study of such regions should provide important information 
on the forming nuclei. In particular we will show that not only the precursor regions act as 
seed for crystal growth, but they also determine which polymorph will be nucleated from 
them. To do so we use the order parameters and Wq, which are very useful in the 
detection of polymorphs. We report their exact definition in the Methods section, and just 
report here their basic properties. is a good order parameter to distinguish between 
bcc crystals and close-packed crystals (hep and/or fee), since it is positive in the former 
whereas negative for the latter. W4 is instead good to distinguish between fee crystals (for 
which it has negative values) and hep crystals (for which it has positive values). Figure |3^ 
shows the probability distribution for the order parameter W4, in liquid regions having Qq 
higher than a fixed threshold, Q^q"^. The PV4 distribution was obtained by considering only 
liquid particles (crystal particles are not included in the histogram) in the metastable state 
(before the critical nucleus is formed), and the Q^^"^ threshold values are always within the 
liquid distribution. While the metastable liquid has on average a symmetrical distribution 
around W4 = 0, Fig. [3^ reveals that the high regions have a predominant contribution 
from negative values of W4, which correspond to the fee symmetry. Similar histograms are 
obtained if instead of thresholds one uses small Qq intervals centered at progressively high 
values of Qe (always within the liquid distribution). 

Since we have shown that crystals form from particles of high Qg? the following scenario 
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emerges for the nucleation of hard-sphere crystals: the supercooled liquid develops regions 
of high orientational order (Fig. Wp), whose symmetry favours the nucleation of the fee phase 
(Fig. [3^). Figure [3)3 plots the probability distribution for the order parameter Wq, showing 
that indeed the regions of high Qq display no preference for the bcc symmetry (characterized 
by Wq > 0). Figure |3t displays the radial distribution function, g{r), for the same high Qq 
regions. Notably, higher Qq regions show an enhancement of the shoulder in the second peak 
of the pair distribution function, which is known to be a structural precursor to the freezing 
transition 5l|. The fact that regions of high Qq are more prone to crystallization can also be 



seen in Fig. |311, where the two-body excess entropy |52|, |53|, S2, is plotted for different values 
of the threshold Qf^'^- It is known that the two-body excess entropy forms the dominant 
contribution to the excess entropy, of the order of 85 — 90% in simple monoatomic liquids. 
Its value is S2 = —6.8 for the metastable liquid, and S2 = —10 for the bulk crystal. The 
inset shows that the S2 value indeed rapidly decreases for increasing values of the threshold 
Q^f^'^. Moreover, the dashed and dotted-dashed lines display the values of S2 calculated for 
particles having < (fcc-like) and W4 > (hcp-like) respectively, demonstrating that 
there is a large difference in the configurational entropy (at the two-particle level) between 
particles having fee and hep symmetry, the former ones being strongly favoured towards 
crystallization (the difference between the S2 value of hep and fcc-like particles is of the 
order of 1%). This implies that although fee and hep have the same free energy in bulk, 
small clusters of fee symmetry have a lower free energy (lower configurational, but higher 
correlational entropy) than those of hep symmetry. 

Discussions 

In this Article we have studied the process of crystallization from the perspective of both 
local translational and bond orientational order. Crystallization has so far been described 
by translational ordering of the density field. However, our study clearly indicates that 
symmetry selection due to packing constraint or directional bonding (like in water, see 
Supplementary Information), which is represented by bond orientational order, plays a key 
role in the crystallization process. It is bond orientational order and neither density nor 
translational order that triggers crystal nucleation. Structuring before densification was 



also reported in Refs. 18|, |25|, ISJ] for molecular liquids, in which the range of the interaction 
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is longer than the size of the constitutive particles. For these liquids, the DFT approach 
showed that the density of the critical nucleus might deviate significantly from the bulk 
phases, and that the compactification of the nuclei has to be accompanied by an increase of 
their local structure (i.e. lattice periodicity). This nucleation mechanism is often ascribed 
to the low compressibility of the liquid, which favours structuring prior to densification 
upon nucleation js]. Contrary to this scenario, our results suggest that this behaviour is a 
consequence of the weak coupling between density fluctuations and bond-order fluctuations, 
with the latter driving the crystallization process. In other words, it is due to the fact 
that nuclei form in precursors of high orientational order, where small displacements can 
considerably increase the order with very little density change (see the pictures in Fig. 
The increase in the structural order parameter is thus inherited from a well defined region 
of the metastable liquid phase space, characterized by having high orientational order. Note 
that high density regions are not necessarily characterized by high orientational order, and 
thus they alone cannot trigger the nucleation process. The existence of these regions of high 
density and low orientational order suggests that it is not the low compressibility of the 
liquid which is responsible for the coupling between density and orientational order. 

Moreover we found that regions of high bond orientational order within the metastable 
liquid not only act as crystal precursors but can also determine which particular crystal 
polymorph will nucleate from them, even when precritical nuclei (which naturally populate 
the metastable phase) are disregarded from the analysis. Since the large population of sub- 
critical embryos belongs to the same metastable free energy basin of the liquid (as can be 
seen for example in Fig. [2b), it is natural to expect that the emergence of polymorphism will 
be a continuous process starting in the liquid phase. Polymorphism develops together with 
bond orientational order, highlighting the role of precursor regions in the polymorph selection 
process. A liquid has orientationally ordered precursor regions which can exist in a variety of 
crystal symmetries, according to some high-dimensional probability distribution. For hard 
spheres, the projections of this probability distribution along some reaction coordinates are 
reported in Fig. |3l and show that precursor regions with the fcc-symmetry are more abundant 
than hep-symmetry regions. So if a nucleation event occurs in any of these regions, the 
crystal environment will reflect the symmetry of the precursor, or the symmetry favoured in 



a liquid. For hard spheres, the preference towards 



experiments j2|, |29|, l34j-l36j and simulations |37H39j , and can be explained classically neither 



cc was pointed out in earlier studies, both 
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by Ostwald's step rule nor by the Alexander- McTague scenario 55| . While correctly pointing 
to the relevance of metastable states, Ostwald's rule cannot be literally applied to predict the 
outcome of a nucleation process. Instead, we have shown that the relative abundance of one 
polymorph over the other depends directly on the liquid-state precursor's composition. This 
may be related to the scenario proposed by Stranski and Totomanow 56], where the embryos 
that form most readily are those with the lowest free energy barrier to nucleation. Our results 
suggest that the physical mechanism behind this rule is a matching of bond orientational 
symmetry between precursor regions and crystals, which leads to the reduction of the free 
barrier for nucleation (the interfacial energy). To give a more quantitative account of this 
scenario, we also calculated the pair correlation entropy of precursor regions (Fig. [3]), showing 
indeed an imbalance between the different crystal symmetries. We confirm that this scenario 
of crystallization and the resulting selection mechanism of polymorphs are also valid for soft 
spheres (the Gaussian Core model) [STj and water (see Supplementary Information). These 
results could help clarifying the mechanism behind the interplay between crystallization and 
liquid polymorphism which was recently found for both water 58| and silicon 59 |. 

Our two-dimensional analysis also unveiled a density range of stability of the crystals 
which continuously form in the metastable liquid. This range of stability is limited at low 
densities by the usual disordered liquid configurations, and at high densities by fivefold 
arrangements of particles. This result, obtained from purely static arguments, provides a 
thermodynamic justification of Frank's hypothesis that icosahedral clusters of particles 
act as inhibitors of crystallization. This finding may enhance our understanding of the 
nature of a supercooled metastable liquid state and crystallization, possibly shedding light 
on the interplay between crystallization and vitrification {goI . Liquid and crystal often have 
very different densities, due to the translational order of the latter. However, our study 
reveals that bond orientational order is the first step in the pathway from the liquid to 
the crystal state, and a disconnection of this link by competing orientational orderings or 
random disorder may be responsible for the avoidance of crystallization, i.e. vitrification 



(see Fig. EJi and Refs. [231, |26|). 
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Methods 



We study the crystallization process in a system of = 4000 mono disperse hard spheres 
of diameter a by means of isothermal-isobaric (NPT) Monte Carlo simulations. Lengths are 
given in units of the particle diameter a and pressure in units of ksT/a^, where ksT = 1. 
We place the spheres randomly in a simulation box at packing fraction rj = 0.5352 and 
equilibrate the system at reduced pressure /3pa^ = 17.0. At this pressure the liquid is 
metastable with respect to crystallization, with a difference in chemical potential between 



the liquid and solid state of /3|A/i| = 0.54. As shown in |39| the free energy barrier between 
the metastable liquid phase and the crystal phase is PAF ^ 18, and the size of the critical 
nucleus is ^ 80. Under this conditions crystallization is a rare event, for which not only long 
trajectories can be obtained for the supercooled liquid, but also enough nucleation events 
can be observed spontaneously. 

To identify crystal particles we use the local bond-order analysis introduced by Steinhardt 
et al. 3l|. One first introduces a (2/ + 1) dimensional complex vector (q;), which is defined 
for each particle i as qim{i) = J2f=i^ ^«m(fij), where I is a free integer parameter, m is 
an integer that runs from m = —I to m = I, Yim. are the spherical harmonics, fij is the vector 
from particle i to particle j, and the sum goes over all neighbouring particles Nb{i) of particle 
i. Since for hard spheres it is known that the stable crystals are the close packed structures 
we can impose N^^i) = 12, i.e. we consider only the 12 nearest neighbours (a procedure 
which is density independent and greatly improves the statistics). /^From the vectors 
one can construct different invariants, and our bond orientational order parameter is one of 



them, specifically qe{i) = \J ^T^Yl^m=-& k6m(0P/(2^ + 1)- The vectors q^ have been proven 
to be useful also to identify crystal particles within the liquid. This procedure, first applied 



to study nucleation by Frenkel and co-workers [28|], consists of comparing the orientational 
environments of two neighbouring particles via a scalar product ■ qeOO/lqeO)!- 

If the scalar product between two neighbours exceeds 0.7 then the two particles are deemed 
connected. We then identify particle i as being within a crystal if it is connected with at 
least 7 neighbours, and otherwise within a liquid. The structural order parameter, S{i), of 
a particle i (which we employed in Fig, [T]3 and Fig. [2}d) simply expresses the number of 
connected neighbours in a continuous way, i.e. S-, = V^'^*) \^'^'f'}'\?''^^t\\ ■ 

^ ^' * ^J=l |qc(«)l|q6U)| 

To distinguish between the different crystal polymorphs we employ the spatially aver- 
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aged local bond order parameters introduced in Ref. 50j. We first define the quantities 



(jimii) = jfji^^k=o qirn{k). Giveu the previous definition, one can construct the rotation- 
ally invariant quantities 



Qi{i) = v/47r/(2/ + l)|qa 

and 



I I I 



I,m2,m3=0 m2 msj l^li^W 

where the term in parentheses is the Wigner 3 — j symbol (which is different from zero only 
when mi + m2 + = 0). 

More details about these analyses are given in Supplementary Information. 
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FIG. 1: Composition and radial profiles for crystalline nuclei averaged over many 
independent trajectories at I3pa^ = 17. a, Relation between cluster size and polymorphs. 
Average number of particles for bcc (circles), hep (diamonds) and fee (squares) polymorphs as a 
function of the total crystal size (n). The dashed line grows as the volume, ~ n, whereas the 
dashed-dotted line grows as the surface, ~ n^^^. The vertical dashed line indicates the critical 
nucleus size ric, which separates the nucleation regime (the tint blue colour region) and the growth 
regime. The inset shows the average density of particles belonging to the different polymorphs, 
and the continuous line the average density of the liquid phase. Also shown are two examples 
of snapshots of crystal nuclei from the computer simulations, at sizes n = 40 (left) and n = 220 
(right). The particles are coloured according to the following code: fee (red), hep (green), and bcc 
(blue), b, Average density profiles as a function of the distance r from the centre of mass of the 
nucleus. Lines are density profiles for nuclei of sizes between n = 5 and n = 205 (plotted every 
An = 20 with the order given by the arrow); each density profile is averaged over nuclei of sizes 
n ± 5. Crystals are nucleated at conditions very far from the bulk value, indicated by the dashed 
horizontal line. The inset shows the density difference Ap between the centre of the nucleus and 
the liquid density, as a function of the normalized nucleus size (-R/-Rcriticai)- c, Comparison between 
the density profile {p{r) black line) and the structural order parameter profile (S'(r) red line) for 
the critical nucleus (size n = 80). Both profiles are normalized to be unity in the fee crystal state, 
and zero in the liquid phase. The inset shows the {S,p)-map for nuclei of different sizes (the same 
as in the panel b). The continuous line is the classical behaviour, while simulation results show 
that nuclei form in ordered precursors, and not locally denser precursors. 
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FIG. 2: Roles of density and bond orientational order in crystal nucleation. a, Rela- 
tion between density (p) and bond orientational order (ge) in the metastable liquid for different 
pressures, (3pa^ = 11, 14, 17. Dashed lines are averages over particles in liquid-like environments, 
whereas full lines-|-symbols are averages over those in crystal-like environments. For each pressure, 
the stable phase is given by the lowest line. For (Spa"^ = 11) which is just below the melting pres- 
sure I3pa^ = 11.54, the liquid line is always stable against the crystal line. But as the pressure is 
increased the crystal line crosses the liquid line to become the stable phase. The transition from 
liquid-like to crystal-like happens at constant density, and can be rationalized by the small cage 
rearrangements (as seen in the snapshots) which are sufficient to promote the transition with very 
little density change. At higher densities a second crossover occurs, and the liquid branch becomes 
stable again against the crystal-like branch, b, Probability density for the structural order param- 
eter S in the {Qe,p) plane. The number of connected neighbours grows continuously from to 12 
from the liquid to the crystal phase. Contour lines are almost parallel to the p axis signalling that 
crystallization is promoted mostly by bond orientational order. Regions of high p contain particles 
in a range of environments from liquid-like to crystal-like, which means that density fluctuations 
alone are not sufficient to promote crystallization. 
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FIG. 3: Mechanism of polymorph selection, a, Order parameter W4 for liquid particles having 
Qe > 0.25,0.26,0.27,0.28,0.29,0.30,0.31,0.32 (the order is given by the arrow). The dashed line 
is the probability distribution for crystalline particles in the same system. As Qq increases, the 
regions of high structural order in the liquid are characterized by a growing population of fcc-like 
clusters, b, Order parameter Wq for liquid particles having Qq > 0.27,0.28,0.29,0.30,0.31,0.32. 
As Qq increases, the distributions move to lower and negative values of Wq, thus showing no 
preference for the bcc symmetry (Wq > 0). Ordering seen in the pair correlation function, c, Pair 
distribution function, g{r), for liquid particles having Qq > 0.25, 0.28, 0.30, 0.32 (the order is given 
by the arrows). The y axis has been split to display the first maximum of g{r) (the corresponding 
X scale is on the top axis). Regions of high Qq clearly show an enhanced shoulder in the second 
peak of the pair distribution function, which is a precursor to crystallization, d, Two-body excess 
entropy S2 (continuous line) , calculated for liquid particles with Qq > Qf^^ ; the dashed and dotted- 
dashed lines are instead calculated for liquid particles having W4 < and W4 > respectively, 
fcc-like particles {W4 < 0) in regions of high Qq are thus favoured for crystallization over hcp-like 
particle (W4 > 0). 
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Supplementary Information is organized as follows. First we discuss in more details the 
methods employed in our study, with emphasis on the crystal identification protocol. We 
then show that the microscopic mechanism of crystallization discussed for HS in the main 
text, also applies to other relevant classes of systems. In particular we examine: i) the 
Gaussian Core model (GCM), a model for self-avoiding polymers in an athermal solvent, 
which belongs to the class of ultrasoft potentials ii) the Mw model, a model for water, 
which belongs to the class of tetrahedrally coordinated liquids 



Crystal identification 

The definitions of the order parameters employed in our study are reported in the Methods 
section. Fig. [S4] shows the distribution of these order parameters for both the supercooled 
liquid and the bulk fee, hep and bcc crystals at reduced pressure f3pa^ = 17. 

The Qa-Qg map (Fig. [S4k ) shows that crystal structures are always located at higher Qq 
than the liquid state. But the map is not effective for distinguishing between the different 
polymorphs due to the large overlap between the bcc and the hep structures. Another warn- 
ing concerning the use of this map for crystal identification regards its pressure dependence. 
As we saw in the main text, crystal particles form at conditions very far from bulk values 
and in particular the density of formation of the smallest nuclei is much lower than the final 
bulk density. By computing the (^4 — Qe map at different average densities for the bulk 
crystals one sees that these maps significantly shift to lower Qq as density is decreased. We 
thus conclude that Q4 — Qq maps cannot be reliably used for crystal identification. Instead, 
to identify the crystal polymorphs we take advantage of the different symmetries that the 
crystals have on the Wq and W4 axis. The bcc structure is in fact characterized by a positive 
Wq distribution (Fig. [S4b ) whereas hep and fee both have negative Wq but differ respectively 
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for their positive and negative values of W4 (Fig. [S4H). We have checked that these symme- 
tries are left unchanged if computing the maps at different average densities (or pressures). 
We finally adopt the following criteria for crystal classification: First crystal particles are 
identified as described in Methods and SI. Then, we identify i) bcc particles as all crystal 
particles with Wq > 0; ii) hep particles as all crystal particles with Wq < and W4 > 0; iii) 
fee particles as all crystal particles with Wq < and W4 < 0. 

Now we justify the claim made in the article that high particles correspond to icosa- 
hedra. To spot icosahedra we follow the definitions in Ref. |3|, which we briefly summarize 
here. Icosahedra can be identified by thresholding the value of the following order parameter 
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Fig. S4: Order parameter maps for the thermal crystals and the supercooled state. 

a, Q^-Qe plane, c, Q^-Wq plane, d, Q4-W4 plane. In b the probability distribution for the 
supercooled state is superimposed on the Qi-Qe map for the perfect crystals (fee, hep and bcc). 
The maps show that crystals have higher values of Qq than the supercooled liquid. We can also see 
that the supercooled liquid prior to crystal nucleation consists of its major liquid portion (the dark 
gray region; Qq < 0.35) and minor solid portion (the tint gray region; Qq > 0.35). The polymorphs 
can be identified by exploiting their symmetries along diff'erent axes: bcc crystals have Wq > 0, 
whereas hep and fee crystals are both characterized by Wq < and have respectively W4 > and 
W4 < 0. 
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Fig. S5: Probability density map in the qe—w^ plane for the metastable fluid at I3pa^ = 17. 

The map clearly shows that high qq particles are characterized by a low value of wq, and are thus 
particles in icosahedral environment. 
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Fig. S6: Free energy barrier for the system at f3pa^ = 13, 15, 17, obtained from Umbrella 
Sampling simulations. 



where the term in parentheses is the Wigner 3 — j symbol (which is different from zero only 
when mi + m2 + = 0), and qim are the Steinhardt bond orientational order parameters 
defined in the Methods section. Following Ref. js] icosahedral particles can be identified 
as particles having wq > 0.023. As shown in Fig. [S5| particles with high gg all lay within 
the region which identifies icosahedral environments. We have thus shown that the high qq 
particles, which form the final liquid stable branch at high densities, are in fact icosahedral 
particles. 

With the criteria for identifying crystal particles it is possible to obtain the free energy 
barrier and the critical cluster size from Umbrella Sampling simulations, where a biasing 
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potential is added to the system Hamiltonian to sample crystalline clusters of large sizes. 
The details of the implementation can be found in Ref. j4|. Fig. [S6] plots the free energy 
barrier /3AF as a function of cluster size. The free energy barrier between the metastable 
liquid phase and the crystal phase at /3pa^ = 17 is /3AF ^ 18, and the size of the critical 
nucleus is ^ 80. These results are in good agreement with the ones in Ref. [s]. In this 
condition crystallization is a rare event, for which not only long trajectories can be obtained 
for the supercooled melt, but also enough nucleation events can be observed spontaneously. 

The composition of nuclei obtained from the spontaneously nucleating trajectories were 
compared with the ones obtained in equilibrium from the Umbrella Sampling configura- 
tions for nuclei of size up to 250 particles. No difference in the average composition of 
the nuclei was found between the Umbrella Sampling configurations and the configurations 
obtained from the Monte Carlo trajectories. This proves that the small clusters are in 
quasi-equilibrium, due to the presence of a free energy barrier. 

Nuclei composition was calculated also for pressures I3pa^ = 13, 15 with Umbrella Sam- 
pling configurations and no sensible change in the polymorph composition was found with 
respect to the reported pressure (3pa^ = 17. 



Gaussian Core model 

The Gaussian Core model (GCM) describes the effective potential between the centres 



of mass of polymers dispersed in a good so 
components, first introduced by Stillinger 



vent. It consists of pairwise sum of Gaussian 
1 

1|. The GCM belongs to the class of ultra- 



soft potentials, for which there is no divergence at contact. Unlike HS, soft particles can 
crystallize in open structures, such as the bcc crystal. We have recently published a detailed 
account on the nucleation pathway in the GCM jo], and thus we report here the results 
relevant for our new analysis. In the following we study the nucleation in the GCM for 
P = 0.05 and T = 0.0052, where the units of length and energy are given by the standard 
deviation and amplitude of the Gaussian potential (as usual in the literature (7|). According 
to the phase diagram calculated in Ref. |8|, the chosen state point has the bcc as the stable 
bulk crystal. We follow the crystallization of 200 isobaric Monte Carlo trajectories starting 
from a metastable fluid phase. Crystal particles are identified with the following set of 
parameters (defined in the previous section), Nc = 9 and gthr = 0.6, and the different 
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Fig. S7: Polymorph selection in the GCM. a, Fraction of particles (/) in a given crystalline 
state as a function of the total crystal size n for P = 0.01. The vertical dashed line denotes the 
size of the critical nucleus ric- b, Density probability distribution. The continuous and dashed line 
are the density histogram for solid and liquid particles respectively. The dots represent the density 
histogram for liquid particles fulfilling the condition > 0.3. c, W4, probability distribution for 
liquid particles having < and Qe higher than a fixed threshold Q^^^- The threshold values 
plotted are Qq > 0.25, 0.27, 0.29, 0.30, 0.31, 0.32. d, We probability distribution for liquid particles 
having Qq > 0.3. 



polymorphs distinguished with the same criteria as HS. 

Fig. [S7k shows the composition of crystalline nuclei as a function of the nucleus size, for 
the hep, fee and bcc polymorphs. The bcc phase is the dominant phase and its fraction 
increases as the crystal nucleus gets bigger. Particles in the fee phase account for ~ 30% 
of the solid particles in the small nuclei, and this fraction decreases as the nuclei become 
bigger. The hep phase accounts only for ~ 20% of solid particles in small nuclei, with 
this fraction steadily decreasing as the nuclei become bigger. The vertical dashed line in 
Fig. IS7k indicates the size of the critical nucleus {ric) obtained from the mean-first passage 
time analysis (see Ref. jsl). The composition of the nucleus for n < is approximately 
constant, whereas, for n > 7ic, the fraction of the bcc polymorph increases at the expenses 
of both the fee and hep phases. 

Fig. [STb shows the density histogram for liquid (dashed line) and solid particles (con- 
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tinuous line). Circles in Fig. [S7b display the density histogram for liquid particles having 
a value of Qq higher than 0.3, showing that it coincides with the density histogram of the 
solid particles. Also for the GCM, high Qq regions are characterized by the same density 
fluctuations as the crystalline particles. 

From Fig. [STk we see that for n < half of the crystalline particles are in the bcc phase, 
whereas the other half are in the fee or hep phase. This is consistent with the Wq map shown 
in Fig. [S7b . which shows an almost symmetrical Wq distribution of liquid particles having 
high Qq. So, unlike HS, the high Qq regions in the liquid phase have a symmetry which 
favours also the bcc phase, promoting its nucleation. Fig. [S7k also shows a fraction of fee 
particles twice the fraction of hep particles. Again this is predicted from the W4 probability 
distribution function in the liquid phase, depicted in Fig. \S7\ i. As in HS, regions of high 
orientational order show a preference for the fee symmetry. 



Mw model for water 

The monoatomic model of water (Mw) is essentially a reparametrization of the Stillinger- 
Weber potential to account for the structural and thermodynamic properties of water 
The model has been very successful in describing the supercooled behaviour of water and, 
unlike all-atom models, it crystallizes relatively easily. Because of its distinctive physical 
properties and its paramount importance, water is a very good test for our microscopic 
description of crystallization. Unlike both HS and GCM, 

• the density of the solid phase is lower than the liquid phase. We should then ex- 
pect an anticorrelation between bond orientational order and density, i.e. the density 
decreasing (instead of increasing, as in HS) with an increase of bond orientational 
order; 

• the solid phases of water at ambient pressure are the hexagonal ice (Ih) and the cubic 
ice (Ic), providing a possibility to test our polymorph selection criteria for these target 
crystals. 



To detect crystal particles we make use of the CHILL algorithm [10|, which is a natural 
extension of the methods described previously for HS to account for the tetrahedral arrange- 
ment of particles. For particle i, we define the four closest neighbours (ji ■ ■ ■ j^) and for each 
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Fig. S8: Polymorph selection in the Mw water, a, Density probability distribution for 
liquid particles having Qq > Qf\ with Qf' = 0.09, 0.10, 0.11, 0.12, 0.13 and the order given by the 
arrow, b, W4 probability distribution for the same set of particles as in panel a. 

bond we compute the scalar product S3 = ^^{i) ■ qsO)- Staggered bonds are characterized 
by S3 < —0.8 and ecHpsed bonds by —0.3 < S3 < 0.1. The Ih crystal is characterized by 
having one eclipsed bond and three staggered bond, whereas Ic crystal has all four bonds 
in a staggered co nfig uration. Also defective crystal configurations are detected by following 
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10|. The definition of these defects varies somehow in the literature, as in 
12| . but we tested that our results are independent of the details of this choice. 
Testing the symmetry of liquid particles is more subtle. In order to find an order parameter 
which can distinguish between the Ic and Ih symmetry, we need to take into account the 
second coordination shell, which comprises 16 molecules. In this case we have found that 
the W4 provides a very good way to distinguish between Ih and Ic, with Ih having W4 > 
and Ic W4 < 0. 

Simulations at both T = 180 K jisl . [l^ and at T = 220 K 12 1 have shown that ice 
spontaneously nucleates preferentially in the Ic form. We run simulations at ambient pres- 
sure and at an intermediate temperature, T = 206 K, and confirm the preference for Ic 
nucleation over Ih nucleation. We then test the symmetry of regions of high bond orienta- 
tional order in the liquid phase (not considering crystalline particles). Fig. [S8k shows the 
density distribution for liquid particles of high bond orientational order, having Qq > Q^f^^, 
confirming that indeed regions of high orientational order are anticorrelated with density. 
Fig. [SSb finally shows the probability distribution on the W4 axis for liquid particles with 
Qq > Qf"^. The distribution becomes more and more peaked towards negative values of 1^4, 
which correspond to the Ic symmetry, as bond orientational order increases in the liquid. 
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This explains why Ic crystals are the most abundant polymorph at this state point. 

We have thus shown that in a large class of potentials (hard, ultrasoft and tetrahedral) 
nucleation occurs always in regions of high bond orientational order, and that these regions 
share the same symmetry of the nucleating solid phase. This suggests the universality of 
our scenario. 
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